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Ensembles of isotropic random matrices are defined by the invariance of the probability measure 
under the left (and right) multiplication by an arbitrary unitary matrix. We show that the multipli- 
cation of large isotropic random matrices is spectrally commutative and self-averaging in the limit 
of infinite matrix size TV — » oo. The notion of spectral commutativity means that the eigenvalue 
density of a product ABC ... of such matrices is independent of the order of matrix multiplication, 
for example the matrix ABCD has the same eigenvalue density as ADCB. In turn, the notion of 
self-averaging means that the product of n independent but identically distributed random matrices, 
which we symbolically denote by AAA. . ., has the same eigenvalue density as the corresponding 
power A n of a single matrix drawn from the underlying matrix ensemble. For example, the eigen- 
value density of ABCCABC is the same as of A 2 B 2 C i . We also discuss the singular behavior of the 
eigenvalue and singular value densities of isotropic matrices and their products for small eigenvalues 
A — > 0. We show that the singularities at the origin of the eigenvalue density and of the singular 
value density are in one-to-one correspondence in the limit TV —¥ oo: the eigenvalue density of an 
isotropic random matrix has a power law singularity at the origin ~ A|~ s with a power s 6 (0, 2) 
when and only when the density of its singular values has a power law singularity ~ A~ CT with a 
power a = s/(4 — s). These results are obtained analytically in the limit TV — > oo. We supplement 
these results with numerical simulations for large but finite TV and discuss finite size effects for the 
most common ensembles of isotropic random matrices. 

PACS numbers: 02.50.Cw (Probability theory), 02.70.Uu (Applications of Monte Carlo methods), 05.40.Ca 
(Noise) 
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I. INTRODUCTION 

Ensembles of Hermitian random matrices with invariant measures have been thoroughly studied in the literature [l|- 
13] • Much less known are non-Hermitian random matrices [j| . In this paper we discuss a class of isotropic non- Hermitian 
matrices that represent a natural extension of the class of invariant Hermitian matrices to the non-Hermitian case: 
the probability measure of an isotropic random matrix ensemble is invariant under the left (and right) multiplication 
by an arbitrary unitary matrix. Here we are interested in properties of the limiting eigenvalue densities of products of 
isotropic matrices in the limit of infinite matrix size TV — ¥ oo. These properties can be deduced from the correspondence 
between large random matrices and free random variables 6], and most of them follow from the Haagerup-Larsen 
theorem Q, that was formulated in the framework of free probability. This theorem gives a very useful relation 
between the eigenvalue density of an isotropic matrix A and the density of its invariant Hermitian partner AA^ . We 
exploit this relation to discuss the spectral commutativity and the self-averaging of the product of isotropic random 
matrices in the large TV limit. The product of identically distributed independent matrices has the same eigenvalue 
density as the corresponding power of a single random matrix [8]. This is an exceptional property which has no 
counterpart in classical probability theory. 

The paper is organized as follows. In Section |TT] we recall the definition of isotropic random matrices and the 
Haagerup-Larsen theorem [7fl. In Section [IIII we discuss products of isotropic matrices and the spectral commutativity 
of the multiplication of infinitely large matrices from this class. In section HVl we illustrate how to use the Haagerup- 
Larsen relation to calculate the eigenvalue density for a few isotropic matrices in the large TV-limit. In Section [V] 
we analyze the correspondence between the singularities of the eigenvalue density and the singular value density. In 
Section PVTl we discuss finite size effects for three generic ensembles of isotropic matrices. In Section fVIII we compare 
eigenvalue densities for products of finite matrices, obtained by Monte-Carlo simulations, with the limiting densities 
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calculated analytically using the method discussed earlier in Section Mil In Section IVIIII we shortly summarize the 
paper. 

II. ISOTROPIC RANDOM MATRICES 

Before we discuss isotropic matrices, let us recall invariant Hermitian random matrices. A random Hermitian 
matrix h is called invariant if its probability measure is invariant under the transformation h — » U~ 1 hU where U is an 
arbitrary unitary matrix. The notion of random matrix is analogous to random variable, so one has to remember that 
the term "random matrix" does not refer to a single instance but to an ensemble of matrices with a given probability 
measure. Thus, invariance of h stands for the invariance of the probability measure or, in other words, that the 
random matrices h and U~ 1 hU have the same probability measures. Randomness of an invariant Hermitian matrix 
is entirely encoded in its eigenvalue distribution in contrast to non-invariant random matrices. However, with any 
Hermitian non-invariant random matrix h n one can associate a unique invariant random matrix h that has exactly 
the same eigenvalue distribution as h n : h = uh n u~ l where u is drawn according to the uniform (Haar) measure from 
the unitary group U{N). 

A random matrix H is called isotropic [8[ if it can be decomposed into a product of an invariant Hermitian positive 
semi-definite random matrix h and a Haar unitary matrix u: H = hu. Clearly h 2 ~ HH^ . Statistical properties of H 
are inherited from its invariant Hermitian partner h. Eigenvalues of h correspond to singular values of H. The concept 
of isotropic random matrices is a natural generalization of the concept of rotationally invariant complex variables that 
can be written as z — re 1 ^, where r is a non- negative real random variable and <j> is a random variable uniformly 
distributed on the interval [0, 2tt). 

One should note that isotropic matrices can be constructed also from non-invariant random Hermitian matrices as 
H = u\h n U2, where iti,U2 are independent unitary Haar matrices and h n is an arbitrary Hermitian positive semi- 
definite random matrix which is not necessarily invariant. In particular h n may be a diagonal random matrix with 
independent identically distributed non-negative real random variables on the diagonal. 

The probability measure of an isotropic random matrix is invariant under the right H — » HU (and left H — » UH) 
multiplication by any unitary matrix U . This property can be used as an alternative definition of isotropic random 
matrices. The eigenvalue distribution puiz) of an isotropic random matrix H is circularly symmetric on the complex 
plane. It depends only on the eigenvalue modulus r = \z\, so it can be written as a function of a single real argument 
Ph(z) = qh (\z\). The radial profile of the eigenvalue distribution of H depends on the eigenvalue distribution ph(X) 
of the matrix h and on the size N of the matrix. In two limiting cases of matrix dimensions, N = 1 and N oo, the 
relation between eigenvalue distribution of the isotropic matrix H and of its Hermitian partner h is explicitly known. 

For N = 1 the random matrix H reduces to a complex random variable, and the matrix h to a real non-negative 
random variable. The relation between the two reads H = he 1 ^ where <j> is a phase uniformly distributed on the interval 
[0, 2tt). The relation between the distributions of the random variables H and h can be conveniently expressed in 
terms of the cumulative density function Fh(x) = J Q ph(\)d\ for h and the radial cumulative distribution function 
for if 

F H {x) = f PH (z)d 2 z = f g H (\z\)d 2 z = f 2-KrQ H {r)dr . (1) 

J\z\<x J\z\<x JO 

The relation between the cumulative distributions for H and h reads Fu(x) = Fh{x). It amounts to 27rrp#(r) = Ph{r)- 
Of course this relation holds only for N = 1. 

Below we use the radial cumulative density function Fh{x) also for N > 1. In this case Fh{x) is defined as the 
probability that a randomly chosen eigenvalue of H lies within distance x from the origin of the complex plane. 

As we mentioned, the relation between the eigenvalue distributions of H and h is also known for N — > oo [7], as in 
this limit random matrices can be mapped onto free random variables 9]: invariant Hermitian matrices to free real 
random variables, and isotropic random matrices to R-diagonal free random variables (Toj . In other words, in this 
case one can use methods of free probability Q to derive the relation. We quote here the result that is known as the 
Haagerup-Larsen theorem Ml. This theorem states that the radial cumulative density function of H can be expressed 
in terms of the S-transforrn0, |Tl| for h 2 = HH^ 

S h 2 (F H (x) - 1) = i . (2) 

The support of the radial profile qh{t) extends from r m i n to r mQX , so that the eigenvalue density of H forms either a 
disk of radius r max if r min = or a ring if r min > 0. The disc (or ring) is centered at the origin of the complex plane. 
The internal radius is r^ in = j pi l (X)X~ 2 dX and the external one = J p/ l (A)A 2 dA. The internal radius is positive 
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r m in > and the support of pn is a ring, if for A — > the eigenvalue density Ph{X) decays to zero faster than the 
first power of A: ph(X) ~ A 1+e , e > 0. Otherwise, the support of the eigenvalue density of H is a disk. To be precise, 
the theorem assumes that the eigenvalue density of h has no isolated point masses (Dirac deltas) in the spectrum. In 
other words, the spectrum of h must be continuous. 

Eq. ([2]) tells us that the cumulative distribution Fh implicitly depends on the eigenvalue density of the matrix h 
through the S-transform S^. So, let us recall what the S-transform is [Tl|. It is defined for an infinitely large 
(N — > oo) invariant Hermitian matrix. Let a be such a matrix. Denoting the limiting eigenvalue density of this matrix 
by p a (A), the S-transform is calculated as follows. First, one calculates the Green's function as the Stieltjes transform 
of the eigenvalue density 

= I p^A 
J i z-X 

G a is a complex function defined outside the support / of the eigenvalue density, which consists of intervals on the 
real axis. The Green's function can be expanded in powers of 1/z, and the coefficients of this expansion are equal to 
the moments of the eigenvalue distribution: 

/*.,*; = / Pa(X)X k d\ , (4) 



i 



for k = 0, 1 . . . and /i a; o = 1- One can alternatively define the moment-generating function: 



M*) = -G a - 1 ■ (5) 

Expanding it in z one obtains an infinite power series 4> a (z) = IZfcLi Pa,kZ k if all moments exist. The S-transform for 
the matrix a is defined as 

z + 1 

Sa{z) = Xa(z) , (6) 

z 

where \a is the inverse of the moment-generating function (p a : 

Xa(<Pa(z)) = <Pa(Xa(z)) = Z . (7) 

The S-transform of the product of independent invariant matrices a, b (free random variables) is equal to the product 
of the corresponding S-transforms (Tlj : 

S ab (z) = S a (z)S b (z) . (8) 

Since the S-transform is a complex-valued function the multiplication on the right hand side of the equation is both 
associative and commutative. This property has deep consequences for products of isotropic matrices in the limit 
N — > oo. In particular, as we discuss in the next section, multiplication of isotropic random matrices is commutative 
in this limit. 

Coming back to Eq. ^ we see that in order to calculate the cumulative distribution of an isotropic matrix H, we 
first have to calculate the S-transform for the matrix h 2 = HH^ . Assume we know the eigenvalue density of h. The 
eigenvalue density of h 2 is 

^(A) = ^=p, l (v / A) . (9) 

Finding the Green's function Gh 2 {z) for h 2 as in Eq. ([3]) and then the S-transform S^(z) as in Eq. ([6]), we obtain an 
explicit equation for Fh (Eg. [5]). Actually, if one eliminates the S-transform from Eq. © one can rewrite the latter 
as 

applying to it Eqs. (|5l6p . This equation has been derived independently in Refs. [TH, [l3j]. This form is however less 
transparent than Eq. which explicitly refers to the S-transform and thus uncovers an important connection to the 
commutative nature of the multiplication of the S-transform (Eq. [5]) that is responsible for the spectral commutativity 
of the isotropic matrices in the large N- limit, as we discuss in the next section. 
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III. PRODUCTS OF ISOTROPIC RANDOM MATRICES 

Consider a product of a finite number, n, of independent isotropic random matrices 

P = H 1 H 2 ...H n (11) 

in the limit TV — > oo. The goal is to calculate the eigenvalue density of P given the eigenvalue densities of the Hi's. 
The product of isotropic random matrices is also isotropic so the eigenvalue density of the product can be determined 
from Eq. ([2]) replacing H by P and h by p. So we have to calculate the S-transform for p 2 = PP^ . As a consequence 
of the multiplication law ( Hp t he S-transform of the matrix for p 2 = PP' with P given by Eq. (jlip can be written as 
a product of S-transforms [HI HH : 

n 

S p2 (z) = l[S h z(z) (12) 

3=1 

for h 2 = HjHj. It is a simple consequence of the associativity and commutativity of the S-transform composition 
rule (0. This means that the equation for the radial cumulative eigenvalue density Fp(x) of the matrix P (fTTj) can 
be expressed in terms of the S-transforms for individual factors in the product (|11[) : 

n i 

Y[S h , (F P (x)-l) = - . (13) 
3=1 

If one writes an analogous equation for the product 

Pit = •ff7r(l)#ir(2) •• • Hir(n)i (14) 

where 7r is an arbitrary permutation of the set {1, . . . , n} one obtains exactly the same equation as in (|13|) 

n 

IJS w (F Pir ( a! )-l) = - J , (15) 

3=1 ' ^ 

but now for Fp^. Therefore the functions Fp and Fp n are identical Fp = Fp r . In other words, the eigenvalue 
distribution of the product of isotropic matrices (|14l) does not depend on the order of matrix multiplication in the 
large TV-limit. 

As we mentioned, the multiplication of infinitely large isotropic matrices has another surprising property. The 
product of n independent identically distributed matrices isotropic that we denote by P = HH . . . H has exactly the 
same probability measure as a power Q = H n of a single matrix [§[ . It was first discovered for the product of Ginibre 
matrices [l|| . For the sake of completness we repeat here the argument given in Ref. [8] . Eq. (TlU)) for the product of 
identically distributed matrices takes the form 

S h2 (F P (x)-l) = -^ , (16) 

which means that Fp(x) — Fh (a; 1 /"). On the other hand, by construction, the eigenvalues of Q are given by 
powers of eigenvalues of H: Xq = A^, so the probability that |Aq| < x is equal to the probability that |A#|™ < x: 
Prob(|A Q | <x) = Prob(|A ff |" < a;). It follows that F Q {x) = F H (x 1 ^) and hence Fp(x) = F Q (x). 

To summarize this section, the product of a finite number of independent infinitely dimensional isotropic random 
matrices is an isotropic random matrix. The eigenvalue distribution of this matrix does not depend on the order 
of multiplication. Moreover, if there are independent identically distributed matrices in the product they can be 
substituted by the corresponding power of a single random matrix from the corresponding matrix ensemble. For 
instance the product ACBACCB has the same limiting eigenvalue density as A 2 B 2 C 3 . It is worth mentioning that 
the effect of self-averaging was also observed for a Wigner type of matrices with independent identically distributed, 
entries belonging to the Gaussian universality class [18j . Intuitively, such Wigner matrices become isotropic in the 
large N limit. 

IV. EXAMPLES OF INFINITELY DIMENSIONAL ISOTROPIC RANDOM MATRICES 



In this section we give a couple of examples of isotropic matrices in the limit N — ¥ oo . We keep the convention that 
isotropic matrices are denoted by capital letters and their Hermitian partners by the corresponding small letters. 
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The first example is a matrix A = au, where u is a Haar unitary random matrix and a is an invariant Hermitian 
random matrix with an eigenvalue density given by the quarter-circle law: 

Pa{\) = , AG [0,2] . (17) 

The matrix a 2 has the eigenvalue density ©: 



^ W = mV' AG [0,4]. (18) 



The Green's function (|3|) of a 2 is 



^ = 5-5^ ■ as) 



Using Eq. (JS|) we find the moment-generating function 



1 



and its inverse function (|7|) 



The S-transform (HI) reads 



6 a2 (z) = — (1- VT^4i) -1 (20) 



S„ a (z) = T ^-. (22) 
1 + z 

Inserting it into Eq. ([5]) we find the radial cumulative distribution for the isotropic matrix A associated with a: 

F A (x) = x 2 (23) 
for < x < 1. The radial profile is qa^) — -^F' A (r) = and hence 

Pa{z) = - , |z| < 1. (24) 

The eigenvalue density p a (x) of the Hermitian matrix a and the corresponding radial profile qa{x) of the eigenvalue 
density of the isotropic matrix A = au are shown in Fig. QJi. Clearly in this case the spectrum of the matrix A is the 
same as for Ginibre matrices [20] |. 

As a second example, we consider a random isotropic matrix B = bu associated with an invariant Hermitian random 
matrix b with the following eigenvalue density 

MA) = -\j{a 2 + -X 2 )(X 2 -a 2 _) , A G [a-,a+] , (25) 

where a± = 1 ± ^fa and < a < 1. For b 2 we have 



(A) = ^ V(A+-A)(A-A_) , A G [A_, A+] (26) 
where A± = (1 ± y/a) 2 . It is straightforward to calculate the Green's function: 

Gp(z) = ^ (z + a- 1 - VC* - A+)(z - A_)) , (27) 



the moment generating function 



a + 1 



.(*) = — (1 - Vd-A + ,)(1-A_,)j - -Z- , (28) 
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its inverse function 



and the .S-transform for the matrix b 2 



*^ = (* + l)(as + l) ' (29) 



Using the Haagerup-Larsen theorem ([2]) we find a very simple equation for the radial cumulative distribution Fb (r) 
of the matrix B 

F B (r) = r 4^4 , re[r ,l], (31) 



where tq = yl — a. It corresponds to a uniform eigenvalue distribution in the ring with, the internal radius r m ^ n — Tq 
and the external one r rnax = 1 with the density 

Pb(z) = 1 2 , |«|e[r ,l] • (32) 

Clearly for a — 1 the previous case is restored. In the remaining part of the paper we choose a = 0.9 for the matrix 
B unless stated otherwise. For this choice r min = v^OTT and ps(z) = for \z\ 6 [\/0T, 1] (see Fig. ((TJd)) 
Next, we define a matrix C = cu with c having the following eigenvalue density 



Pc(A) = 27V^T ' Ae[M ■ (33) 



Repeating all steps as in the previous cases we find 



and the S-transform 



4 ( *{z) = -=J^=-\ , (34) 
y/y/l- 16z + l 



M*) = (— !— - , * ) . (35) 
4z \z + 1 (z + l) 3 ' 



For this S-transform the cumulative distribution function for C ([2]) is given by the solution of the following equation: 

4F%(x) - x 2 F c (x) - x 2 =0 . (36) 

This equation can be solved for Fc(x). The solution has three branches, and one has to choose the branch that gives 
a monotonic function on x € [0, y/2] increasing from Fc(x = 0) = to Fc{x = v2) = 1- The radial profile of the 
eigenvalue distribution is obtained by differentiating the solution qc{t) — ~^F' c (r). It is shown in Fig. ((TJ;). 

We additionally consider an isotropic random matrix D = du constructed from an invariant Hermitian random 
matrix d that has a uniform distribution pd{^) — 1 for A € [0, 1]. The radial profile of the distribution for the matrix 
D is shown in Fig. (JTJl). 

In Section [VIII we employ a finite size version of the aforementioned classes of random matrices to discuss products 
of large but finite isotropic matrices. 



V. SINGULAR VALUES OF ISOTROPIC MATRICES 



As follows from Eq. (|2|), the eigenvalues statistics of infinitely large isotropic matrices is in one-to-one correspondence 
with the statistics of singular values. Indeed, Eq. ([2]) provides a relation between the eigenvalue density of an isotropic 
matrix Fh{x) and the S-transform for the matrix h 2 = HH^ . Obviously, eigenvalues of h are equal to singular values 
of H. From the S-transform one can derive the eigenvalue density of h 2 . For instance, for the Ginibre ensemble the 
eigenvalues are distributed uniformly on the unit disk. The eigenvalue density is p(z) = l/ir for \z\ < 1 and hence 
Fh{x) = x 2 for x £ [0, 1]. Inserting this result into Eq. ^ one finds Sj^^z) = 1/(1 + z). This also means that the 
eigenvalue distribution of h 2 = HH^ is given by the Wishart distribution (JT5J) and the one of h by the quarter-circle 
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s Az) = tt-t^;, ( 38 ) 



law (I17p . This is of course the same calculation as the one presented in the previous section for matrix A, but now it 
has been carried out in the opposite direction. 

An interesting situation is encountered for the matrix P = HH . . . H being a product of n independent infinitely 
large Ginibre matrices H. As discussed in Section QUI F P (x) = F H (x 1/n ) = x 2 / n (see Ref. @]), and hence 

PP (z) = -\-FU\z\) = —M- 2 ^ (37) 
2n\z\ nn 

on the unit disk \z\ < 1. The S-transform for p 2 = PP' is 

1 

(1 +"«)' 

as follows from Eq. ([2]). We can now write an explicit equation for the moment-generating function (J7J 

z(l + p2 (z))" +1 =^{z) , (39) 

and for the Green's function 

z n G r ;+ 1 (z) =zG p2 (z)-l . (40) 

It is worth noting that one can derive an analogous equation also for the product of rectangular Gaussian random 
matrices [HI, HH • The solution of Eq. (|39p can be written as a power series 

with the coefficients fit given by the Fuss-Catalan numbers [l7| • The radius of convergence of the power series <p p 2 (z) 
is equal to 

r = ij m _ L_2_ — (42} 

c fe^oo fi k+1 (n + l) n+1 1 v ; 
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and <f> p 2{z) is singular when z approaches r c . The first moment of this distribution corresponds to the mean square 
of singular values of the matrix P: (Xgy)p = fit = 1. This can be compared with the mean squared absolute value 
of eigenvalues of the matrix P, which is given by the second absolute moment of the distribution (|37]l : (\\ev\ 2 )p = 
/ \z\ 2 pp(z)d 2 z = l/(n + 1). We see that (\Xev\ 2 }p < (A|y)p, as expected on general grounds. 

Let us now discuss the behavior of the singular value density of P. Singular values of P are equal to eigenvalues of 
the Hermitian matrix p, so this density is equal to the eigenvalue density p p (X). First, let us determine the behavior 
of p p 2(X) for A — > 0. As follows from Eq. (|40|) . the Green's function has a singularity for z — > 0: 

G p 2(z) ~ Z -^t . (43) 

From Eq. ©, this means that also the eigenvalue distribution must have the same singularity 

,V(A)~A-^ (44) 

when A — > 0, which corresponds to the following singularity of the eigenvalue density of p: 

Pp (X) ~ A"S£ . (45) 

The powers in Eq. (|37l) and Eq. (j45|) are related to each other. For n = 1 both functions are regular at zero (as they 
are given by the quarter-circle and Ginibre laws, Eqs. (ITT)) and (f2"4"|). respectively). For n — 2 the eigenvalue density 
behaves as \z\~ x while the singular value density as A -1 / 3 , for n = 3 they behave as |z|~ 4 / 3 and A" 1 / 2 , respectively, 
etc. One can repeat the discussion also for products of rectangular matrices [Til ITsT]. 

The function G p 2(z) has a cut along the interval [0, (n + l) n+1 /n n ] on the real axis that corresponds to the support 
of the eigenvalue distribution p^ (A) . The upper end of this interval is equal to l/r c (1421) as follows from the change 
of argument z — > 1/z in G(z) and <fi(z) (see Eq. (|SJ)). The position of the upper end of the interval can also be found 
directly from Eq. (|4T))) as a place were the function G = G p 2{z) is singular. Singularity means that either dG/dz = 
or dz/dG = 0. Writing Eq. (|4T)|) as z n G n+1 = zG — 1 and differentiating both sides with respect to G and setting 
dz/dG — we obtain (n + l)z n G n — z. Solving these two equations we find z = (n + l) n+1 /n n , that corresponds to 
the singularity at the upper end of the support of the density p p 2 (A). One can also determine a closed-form expression 
for p p 2 (A) in terms of special functions ■ 

The relation of the singularities for A —¥ is actually more general. For any infinitely dimensional isotropic matrix 
H whose eigenvalue density has a power-law singularity 

Ph{z) ~ \z\- s (46) 

with < s < 2 we can determine the power of the corresponding singularity of the singular value density. From Eq. 
@ we can derive the behavior of the S-transform for h 2 for z — > as Sh?(z) ~ (1 + z)~ 2 /( 2 ~ s \ which implies (see 
Eqs. H])-®) that the Green's function of h 2 has a sing ularity G h 2(z) - z ~ 2 ^ 4 -^ at the origin. This singularity is 
linked to the singularity of the eigenvalue density /^(A) ~ \- 2 /( 4 - s ) for A — > 0. By changing variables A — > A 2 we 
obtain a singularity of the eigenvalue density of h (which is the density of singular values of H ) : 

p h (X) ~ . (47) 

To summarize, an isotropic matrix whose eigenvalue density has a singularity \z\~ s at the origin of the complex plane, 
has a singularity \~ s /( 4 ~ s "> for A — > in the singular value density. This statement can be inverted. An infinitely 
large isotropic matrix H having a power- law singularity in the density of singular values Ph(X) ~ \~ a with < a < 1 
has a power law singularity in the eigenvalue density with the following power 

p H (z)~\z\-&. (48) 



VI. FINITE SIZE EFFECTS 



So far we have discussed the limiting densities for N — > oo. However, in many practical problems one encounters 
large but finite matrices. The calculation of the eigenvalue density for finite N is much more complicated. Moreover, 
contrary to the limit N — > oo, the results for finite N are not universal and they depend on many details of the 
probability measure. The finite iV-density has been calculated analytically only for a couple of very specific cases in- 
cluding Ginibre matrices 0, H^, [2l| , elliptic matrices [22| , unitary truncated matrices [23[ and products of independent 
Ginibre matrices [26l |. 
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Generally, various classes of isotropic random matrices may have the same limiting eigenvalue density for N — > oo 
but completely different properties for finite N. In this section we discuss the three most common classes of isotropic 
random matrices. The first class is isotropic random matrices defined by the partition function EH 

Z = J DH exp (-JVTr V (HH f )) , (49) 

with a potential V(x) which is an TV-independent polynomial (or power series) in x. The symbol DH denotes the flat 
measure for N x N complex matrices. 

The second type are random matrices constructed as 

H' = hu (50) 

where u is an N x N Haar unitary matrix and h is an N x N Hermitian matrix generated from an invariant ensemble 
defined by the partition function 

Z = J Dhexp(-NTrv(h)) . (51) 

The potential v(x) = v(—x) is an N -independent even polynomial (or power series) in x and Dh is the flat measure 
for N x N Hermitian matrices. 

Finally, the third type are random matrices constructed as 

H" = U\hu2 , (52) 

where u\ and ui are independent matrices uniformly distributed (according to the Haar measure) on the unitary 
group U (N) , and h is an N x N diagonal Hermitian matrix. The entries on the diagonal are independent identically 
distributed real (non- negative) random variables with an iV-independent probability distribution, whose density we 
denote by p(x). As far as the eigenvalue spectrum is concerned random matrices (|52l) have the same eigenvalues as 
random matrices defined as a left multiplication of a diagonal Hermitian matrix by a Haar unitary matrix hu. Such 
matrices were studied in Ref. [24[ where they were called sub-unitary. 

In all three cases the probability measures are invariant under the left (H — > UH) and right (H — > HU) multipli- 
cations by an arbitrary unitary matrix U. A common feature of these ensembles is also that the functions p, v, V, 
which define the probability measures, do not depend on N. Otherwise, the three ensembles are completely different 
and have different eigenvalue statistics. For the matrices of the third type the eigenvalues are independent while for 
the other two they are not 

For the sake of illustration let us discuss three ensembles of isotropic matrices, all having the limiting density for 
N — > oo: ph{z) = 1/tt on the unit disc. The matrix H of the first type (|49p is defined by the partition function 
with the potential V(x) = x. It is a Ginibre matrix [2(j. The matrix H' of the second type is defined by the 
partition function with the potential v(x) — x 2 /2. The matrix H" of the third type (|52p is defined by the probability 
distribution with the probability density function given by the quarter-circle law p(x) — \/A — x 2 /ir for x G [0, 2]. In 
the limit N — > oo the eigenvalue densities of H, H' and H" tend to the same limiting distribution equal to 1/tt on 
the unit disc, but for finite N the eigenvalue densities of H, H' and H" are different. For N = 1 they are 

P(z) = -e- |z|2 , (53) 

for the first type, 

p'(z) = e-l z ' 2 / 2 , (54) 

for the second one, and 



p (z) = ( } 

for the third one, respectively. The three expressions can be easily derived since for N = 1 random matrices reduce 
to scalar random variables. The characteristic l/\z\ behavior for H' and H" is generated by the Jacobian of the 
transformation to polar coordinates of the flat measure d 2 z — rdrdcft: since tf> and r are independent random variables 
we have p(r)dr^dcj) = p(\z\)^-rd 2 z. 
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FIG. 2: (a) Numerical eigenvalue densities (for matrix size N = 25 (green crosses) and 400 (blue circles), made of 10 7 
eigenvalues) for Ginibre matrices described by Eq. ((49}, compared with the corresponding N = 1 (Eq. (|53[) , (red solid line), 
N = 2 (Eq. l|A5f)) (red dashed line), iV = 3 (red dashed-dotted line), iV = 4 (red dashed-triple dotted line) and N — > oo (black 
line) densities, (b) Same plot as (a) for matrices of the type described by Eqs. (|50[) an d (|51[) . (c) Same plot as (a) for matrices 
of the type described by Eq. (|52[) . 

One can also calculate the eigenvalue densities for N — 2,3,4,.... The details of the calculations are given in 
Appendix El Except for the Gaussian ensemble of type one (j49|) . where a closed formula for finite A" density can be 
given in a simple form (|A5|) . the calculations for two remaining types of matrices get more and more tedious and 
cumbersome with increasing N. They become easy again for large N, for which one can anticipate a form of finite size 
corrections to the limiting density. In Fig. [2] we show densities for A~ = 1,2, 3,4, oo which are obtained analytically 
and for A" = 25,400 which are obtained by Monte Carlo simulations. The finite-size corrections close to the edge of 
the support of the limiting density at A = 1 take a form of a sigmoidal function. The extent of the crossover region 
of this function shrinks with N to zero in the limit N —> oo. For the Gaussian ensemble of type P9")) the shape of 
the sigmoidal function is known to be given by a complementary error-function with the range parameter scaling as 
1 / vN [5( . Also for the two remaining types of finite-size ensembles (|50l) and ((52")) one expects the crossover behavior 
to be controlled by an error- function with the same 1 / y/N -scaling [iij . Indeed this is what we observe numerically. 

There is a significant difference between the approach of the finite A^ densities to the limiting density at zero for 
the first type of matrices (|49p for which the approach is uniform and the two remaining types (|50I52I) for which it 
is non-uniform (see Fig. [2]). The non- uniform behavior is a remnant of the l/\z\ singularity in Eqs. (1541) and ((551) . 
When N increases the singularity is pushed towards zero and it eventually disappears in the limit. To compensate for 
the excess of eigenvalues in the region close to the origin the finite A" profiles develop a shallow dip for intermediate 
values of A as can be seen by eye for A* = 25 and A" = 400 in Figs. (5]b andfjje. The effect of the excess of eigenvalues 
at the origin of the complex plane is not present for matrices of type one (|49[) because of the repulsion of eigenvalues 
from the origin plj |. 

VII. PRODUCTS OF FINITE MATRICES 

In this section we study how the commutative and self-averaging properties of finite isotropic matrices set in when 
N increases. To this end we exploit finite size versions of A, B, C, D matrices introduced in Section IIVI As A we 
take Ginibre matrices [2(| which belong to the first class ((49)) discussed in the previous section, while as B,C,D we 
take matrices ((52p constructed from diagonal random matrices by isotropic unitary randomization. First we study 
the size dependence of the eigenvalue densities for products of two matrices similarly as we did for single matrices 
in Section IVI1 The finite A" spectra are generated by Monte-Carlo simulations and they are compared in Fig. [3] to 
the corresponding limiting densities for N oo which were obtained from the S'-transform manipulations and the 
Haagerup-Larsen theorem, as described in Section 1111)) . The size dependence of the finite N densities is very weak in 
the bulk of the distribution as one can see in Fig(3J where the data points for A~ = 25, A" = 100 and A~ = 400 lie 
on top of the limiting curve. A significant dependence on N is observed only in the region close to the edge of the 
distribution. In this region the density takes the form of a sigmoidal function. The range of this function shrinks 
when N increases to eventually restore a sharp threshold at the edge in the limit N oo. This effect is analogous 
to that discussed for products of Gaussian matrices in Refs. [14l - [l6j . where it was conjectured that the sigmoidal 
function in that case was given by the complementary error function. 

In turn, in Fig. 2]we restrict ourselves to matrices with N = 100 but we compare densities for products of isotropic 
random matrices multiplied in different order, as for instance ABDC and ACDB. We see on each plot in Fig. 0] 
that data points representing different order of multiplication lie on top of a master curve within the symbol size. 
This means that already for matrices of size of order A~ = 100 the multiplication of isotropic random matrices can be 




FIG. 3: Monte Carlo eigenvalue densities of products of two isotropic random matrices compared with the theoretical predictions 
(black lines) obtained from the solution of Eq. (|13|l . Different plots refer to the products AB (a), AD (b), and BD (c). In 
all three cases, the results obtained with three different matrix sizes (TV = 25 (blue crosses), 100 (red circles) and 400 (green 
triangles)) are shown. All Monte Carlo densities are made of 10 7 eigenvalues. 




FIG. 4: Monte Carlo eigenvalue densities of products of isotropic random matrices. Different plots refer to the products: (a) 
ABC (red circles) and ACB (blue crosses ), (b) ABAC (red circles), A ABC (blue circles) and A 2 CB ( green triangles), and (c) 
ABCD (red circles), ABDC (blue crosses) and ACDB (green triangles). In all plots black lines represent the limiting result 
(TV — > oo) obtained as solution of Eq. (113[l , All Monte Carlo densities have been produced by numerical diagonalization of 10 
matrices of size TV = 100. 



treated in practical applications as spectrally commutative. 

We also check self-averaging by comparing products containing multiple uses of a single matrix to products of 
independent matrices. For example, we consider a product of the type A 2 BC and AABC, in the first of which A is 
used twice and in the second of which two different A's, representing identically distributed but independent matrices, 
are used. Again the deviations between the resulting spectra are small and practically not detectable in the observed 
resolution (see Fig. U]b). However, one can generally observe in the finite TV eigenvalue spectra that the multiple 
use of a single matrix in a product has a stronger influence on the shape of the spectrum than the change of matrix 
ordering in this product. The effect has been studied quantitatively only for products of Ginibre matrices, for which 
one can analytically derive a finite TV eigenvalue density for the product of n independent matrices AA . . . A and for 
the corresponding power A n of a single matrix (26|. In this case one explicitly sees that the range of the sigmoidal 
function corrections at the edge of the spectrum is of order a/\/~N with a coefficient a that increases as yfn for the 
product of independent matrices and as n for the n-th power. Thus the finite size effects are bigger in the latter case. 



VIII. DISCUSSION 



We have shown that the eigenvalue densities of products of isotropic random matrices do not depend, in the large 
TV limit, on the multiplication ordering. They only depend on the random matrix ensembles from which the matrices 
in the product are generated. We have also extended the result on self-averaging Q by showing that a single matrix 
from an isotropic random matrix ensemble is representative enough to describe multiple independent occurrences 
of matrices from the same ensemble within the same matrix product. We have also derived a relation between the 
exponents which determine the singular behavior at zero of the eigenvalue density and the singular value density of 
infinitely large is otropic random matrices. This result generalizes a previously known relation for products of Gaussian 
matrices [14l4l6l |. 

Isotropic random matrices are a very special class of non-Hermitian random matrix ensembles. Commutative and 
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self-averaging properties of products of such matrices in the large N limit follow from the Haagerup-Larsen theorem 
and the existence of a one-to-one correspondence between invariant large matrices and free random variables, as 
well as from the commutative and associative properties of the multiplication law (jHJ. It would be very interesting to 
work out similar relations for products of generic non-Hermitian matrices. The first step in this direction has been 
done in Ref. [28| where the corresponding multiplication law for a large class of non-Hermitian matrices has been 
derived in the large N limit. The generalization of the R and S transforms leads in that case to quaternion-valued 
functions of quaternion-valued arguments. The corresponding multiplication law is not commutative and thus the 
problem is more complicated. 



Appendix A: Calculations of eigenvalue density for finite N 

In this appendix we discuss how to calculate eigenvalue distributions for finite N isotropic random matrices of the 
three types of ensembles introduced in section (|VI[) . For the sake of illustration we concentrate on matrices that for 
TV — > oo have the limiting density p(z) = 1/tt on the unit. 

For matrices of the type (|49| the corresponding matrix H is an N x N complex matrix. The partition function has 
in this case a linear potential V(x) = x 

Z = [ DHe- NTl HH " . (Al) 



This is the standard Ginibre ensemble [20]. Matrices of this type have been thoroughly studied in the literature. 
Here we recall the main results and refer the reader to Ref. ([Ej]) for details. As it stands, the integrand defining the 
partition function depends on A^ 2 complex degrees of freedom corresponding to the matrix elements. When one is 
interested only in quantities depending on the eigenvalues, one can reduce the complexity of the problem by leaving an 
explicit dependence only on the N complex eigenvalues Zi, i = 1, N of the matrix by integrating out the remaining 
degrees of freedom. This gives, up to a normalization constant: 

dz\ . . .dz 2 N P{z\, . . . ,z N ) , (A2) 
where P(zi, . . . , zn) is the eigenvalue joint probability density function 

atN(N+1)/2 

P(z u ...,z n )= N n e -*Sg-iW I] \*i-«\ 2 - ( A3 ) 

n llfc=l fc - l<i<j<N 

It is normalized in such a way that J dz\ . . . dz 2 N P{z\, . . . , Zjj) = 1. Integrating all but one eigenvalue from the joint 
probability density function one obtains the eigenvalue density of H for finite N 

Pn(z) = J dz%...dz% P{z,z 2 ,...,z N ) . (A4) 

Note that P(z\, . . . , zn) is symmetric with respect to permutations of eigenvalues. Therefore it does not matter which 
N — 1 eigenvalues are integrated out to obtain the density. The result reads 

Mz)=e -^ N fiMr. (A5) 

A — ' n\ 

The radial profiles of these functions for different values of N are plotted in Figj2]a. For large N the shape of the 
radial profile is well approximated by a complementary error function changing from to 1 in a region whose size 
scales like 1/yN. In the limit N — >• oo this region shrinks to one point, so the function becomes a step function. For 
a detailed discussion we again refer the reader to the excellent review [5[ . 

Let us now discuss matrices of the third type (|52[) constructed from diagonal matrices h by unitary randomization 
H = U\hu2 with two independent Haar unitary matrices u\ and u 2 - The ensemble of such matrices has the same 
eigenvalue density as the ensemble of matrices H — hu randomized only on one side. The eigenvalue density for the 
latter one can be calculated as 

PB(z) = o ? y:s^(z-x i )) (A6) 

/f . U 
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where A,'s are eigenvalues of H, and the average is taken over h and u. One can first average over it's and then over 
h's: 

Ph{z) = J dhi . . .dh N P(hi : . . . ,PN)ph u ...,h N ( z ) ■ ( A7 ) 
Here ph!,...,h N (z) is the result of averaging over u for a fixed h = diag(/ii, . . . , /in): 

^...^h/^E^z-A,)) ■ (A8) 

\ i=l In 

The integration measure J dh\ . . . dftjvP(/ii, . . . , ft at) is the probability measure for the diagonal matrix h. For inde- 
pendent and identically distributed h's it factorizes as P(h\, . . . .hpf) = p(hi) . . .p(hpj), where is the probability 
density function for diagonal elements. In our case it is p(x) = \/A — x 2 /tt for x € [0,2], so we have 



M*)= I ^l — f 2 dh N ^ 4 K p ( z ) . (A9) 

The eigenvalue distribution p^, ...,/»# (2) has been explicitly calculated in Ref. [27| . In order to write down the result 
it is convenient to order the /ij's: < fti < h 2 . . . < h,N- With this ordering the density ph 1: ...ji N {z) is given by 



1 N 

Ph 1 ...h N (z) = — J] ^(M*!) , ftfe < M < ftfc+i (A10) 



i=fe+i 



with the F's defined on rings \z\ E (hk, hk+i) for k — 1, . . . , N — 1. Inside the disk or radius /ij < fti) and outside 
the disk of radius {\z\ > /ijv) the density vanishes Ph!...h N {z) — 0. The function i* 1 reads 

/I.2 _ |.|2\Af-2 N - 1 i 

f |,|) = ^ 4i W~ 2(w) 7iv=iv + - 1 - > ( An ) 

lljyi^i V i=0 ( e ) 

and the symbols sf j are defined as symmetric polynomials of order ^ in ft| variables, putting hi — 0. In other words, 

one can write s° = 1, s 1 = $^i=i s2 — Si<j ^i^f e * c -> an< ^ take s (i] = s£ Ui=o- For different orderings of the h's, 
result can be mapped onto the one given above by an appropriate permutation of the indices that organizes the h's 
in increasing order. Therefore it is sufficient to calculate the integral (IA9|) only for h i < hi . . . < hjy since it takes the 
same value for all remaining orderings (permutations): 

p H {z)=N\ f^lJ^h\... f ^ILM-hlp hl _ hN {z) . (A12) 

One can also give an explicit form of the function F for the case when two or more h's have the same values [27| . 
but this case is irrelevant from the point of view of the last integral, since it gives a contribution of measure zero. 

For illustration let us give an explicit form of Ph\,h2{z) for N — 2 that follows from (|All|) . Assuming h\ < h 2 the 
density reads 

hi < \z\ < h 2 / A1 g, 
otherwise . 

One can write explicit expressions also for N = 3, 4, . . . using (|A11|) and integrate them over h's (IA12[) . We have done 
that for N = 2, 3, 4. The results are shown in Fig. [5Jc 

Let us make a general remark. We see that already for finite N the matrix H — hu, where u is a Haar unitary 
matrix on U(N) and h is a constant matrix h = diag(/i l5 . . . , h^f), has many properties that are expected in the large 
N limit. It has a spherically symmetric eigenvalue density on a ring whose radii depend on h's. Actually one can 
show [HJ that Eqs. (IA10IA11[) reproduce the Haaregup-Larsen equation in the large N limit, when the distribution 
of h's becomes a continuous function. 

We can apply a similar strategy to the finite N ensembles of matrices H — hu of the second type (|50"|) where now h 
is a Hermitian matrix from an invariant unitary ensemble. We can use again Eq. (|A7|) but with P(h\, . . . , /ijv) being 
the joint probability function [l[ 

P(h 1 ,...,h N ) = C N e-f^ h l Y[(h k ~h 3 ) 2 , (A14) 

j<k 
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with CV being a normalization constant such that J dhi . . . h]yP(hi, . . . , hjy) — 1. There are two essential differences 
with respect to the previous case: the joint probability P cannot be factorized, and the arguments hi of P may take 
negative values. While doing the integration in (|A7j) over h's, it is convenient to restrict to non-negative semi-axes 
hi > 0, for all i — 1, . . . .N. To this end we introduce a new function Q defined for non- negative h's which is obtained 
from P by integrating out signs of h's: 

Q{hi,...,h N ) = P(sihi,...,s N h N ) . (A15) 

Sl=±l,...Siv=il 

Since the density Phi...h N ( z ) depends only on the absolute values of h's {ph 1 ...h N { z ) = P\h 1 \...\h N \( z )) we have 

/OO P OO /"OO /"OO 

dhi... I dh N P(h 1} . . . ,pN)ph 1 ,...,h N (z) = dhi... dh N Q(hi,...,h N )p hu ,„j lN (z) . (A16) 
-oo J -oo JO JO 

For illustration, let us write it for N = 2 

Q(h 1: h 2 ) = 2C 2 e~ h *- h i ((h 2 - hi) 2 + (h 2 + h,) 2 ) . (A17) 
The integral over all positive h's can be now reduced to an integral over ordered sets h\ < h 2 < . . . < as before, 

/>oo />oo 

Ph(z) = N\ dhi... dh N Q(hi, . . . ,pN)ph 1: ...,h N (z) , (A18) 



with Ph!....,h N {z) given by Eq. (|A10[) . Using this method we have calculated ph(z) for N — 2,3,4. The result is 
presented in Fig. [2]b. 
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